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The firing statistics of Poisson neuron models driven by 
slow stimuli 
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Abstract The coding properties of cells with different 
types of receptive fields have been studied for decades. 
ON-type neurons fire in response to positive fluctu- 
ations of the time-dependent stimulus, whereas OFF 
cells are driven by negative stimulus segments. Bipha- 
sic cells, in turn, are selective to up/down or down/up 
stimulus upstrokes. In this paper, we explore the way 
in which different receptive fields affect the firing statis- 
tics of Poisson neuron models, when driven with slow 
stimuli. We find analytical expressions for the time- 
dependent peri-stimulus time histogram and the inter- 
spike interval distribution in terms of the incoming sig- 
nal. Our results enable us to understand the interplay 
between the intrinsic and extrinsic factors that regulate 
the statistics of spike trains. The former depend on bio- 
physical neural properties, whereas the latter hinge on 
the temporal characteristics of the input signal. 

Keywords Sensory systems • Receptive field • Poisson 
neural models • Slow stimuli • Peri-stimulus time 
histogram • Inter-spike interval distribution 

1 Introduction 

Different sensory neurons have different filtering prop- 
erties. The population distribution of these filtering prop- 
erties has been most extensively studied in peripheral 
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sensory neurons. F or example, in a study carried out by 
(jSegev et al.ll2006l ) , the receptive fields of ganglion cells 
in the salamander retina were characterized by their po- 
larity (ON, OFF, biphasic ON or biphasic OFF), and 
temporal scale (fast, medium, slow). A broad distribu- 
tion of receptive field types was found. However, not all 
possible receptive fields were equally represented in this 
distribution: 93% of the studied cells clustered into a 
few well defined prototypes. If there are discrete classes 
of cells, then presumably the nervous system requires 
specific pre-processing operations on the external sig- 
nals. Such operations are carried out by the particu- 
lar types of receptive fields found in each animal. Here 
we aim at providing a rigorous characterization of the 
type of information encoded by single cells, depend- 
ing on their filtering properties. We work with arbi- 
trary Poisson-like neurons, described by a sequence of 
Volterra filters. In those cases where the neuronal dy- 
namics reduces to a cascade model, the Volterra filters 
may be replaced by a linear receptive field, and a static 
non-linear function governing spike generation. Our aim 
is to understand how the filtering neuronal properties 
(captured by the Volterra filters, or alternatively, the 
receptive field and the non-linearity) determine which 
stimulus features are represented in the peri-stimulus 
time histogram (PSTH) of a neuron. In addition, we 
also study how they affect the inter-spike interval (ISI) 
distribution. We derive analytical expressions for the 
PSTH and the ISI distribution, and we discuss how the 
statistical properties of the responses depend on intrin- 
sic and extrinsic factors, differentiating in particular in- 
trinsic cellular properties from stimulus characteristics. 

We work in the limit of slowly varying input signals. 
Our results are thus only relevant for neurons driven 
by slow stimuli. Examples can be found in hippocam- 
pal interneurons receiving massive theta input from the 
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septum (jToth et al.lll997h or in ol factory neurons mod- 
ulate d by the respiration cycle (IFontanini and Bower 
20061) . An even slower action is exerted by modulatory 
neurotransmitters in the br ain. These s i gnals contro l 
the overall level of arousal ( Sapei 200d Dova 20021 ). 
and they affect the intrinsic neuronal properties of their 
target neurons in time scales of seconds, that is, much 
longer than the time scales governing spike dynamics. 
As a final example, slow adaptation processes can also 
be described in term s of intrinsic transient currents with 
long t ime constants ( ErmentroutlllQQSt Benda and Herd 



l2003t IWilson et al. I l2004h . that either increase (facilita- 
tion) or decrease (adaptation) the excitability of the 



symbols we may replace the functions s{ti) by their se- 
ries expansions around the limit of integration t. Thus, 

for ti G [t- Tjn,t], 



A 

4U) = y -f (ii - ty, where 
dP 



(2) 



A, 



2 Derivation of the PSTH 



and Aq = sit)- If the stimulus varies slowly, then Eq. ^ 
is a rapidly converging series inside [t — Tm, t], so we only 
keep the first few terms in the expansion. The parame- 
ters Xj capture the effect of the stimulus on the PSTH. 
In order to separate this effect from the one produced by 
intrinsic neuronal properties, we define the coefficients 
jj]i,---,3k ^ gunimarizing the neurons's filtering character- 
istics, namely. 



In Poisson neurons, the probability of generating a spike 
at time t conditional to a time dependent stimulus s 
may be expanded in a Volterra series as 



/ + CX3 
dtl hl{t-tl) s{ti) + 
- OO 

\-00 /' + 00 



(1) 



/-\-oo f-\-00 
dti / dt2 h2{t-ti,t-t2) s(ii) s{t2) + --- + 
-OO J — OO 

/ + CXD r-\-oo 
dtl... dtk hk{t-tT_,...,t-tk) S{tl) . . . S{tk) + 

- OO J — OO 



Equation ^ is valid under the assumption that the 
firing probability only depends on the stimulus history, 
and not on the previous activity of the neuron. As a con- 
sequence, the kernels hk only depend on the differences 
t — ti, as opposed to separate dependencies on t and ti. 
Thus, within this picture, intrinsic bursting, refractori- 
ness, adaptation, and all response-induced modulation 
processes are neglected. 

The kernels h^ixi, . . . ,Xk) quantify the efficacy of 
s(— xi), . . . , s(— Xfc) in modulating the spiking proba- 
bility at time t — 0, and ho accounts for spontaneous 
firing. Each hk is measured in units of [time]^^'^+^^ • 
[stimulus] Causality imposes that hk{xi, . . . ,Xk) — 
0, if any Xi < 0. 

Stimuli occurring in the distant past cannot affect 
the firing probability at present times. Hence, we de- 
fine Tm as the memory time of the system, such that 
all kernels hk{xi, . . . , Xk) ~ 0, if any of the Xi > t^. 
Hence, the integration limits (— oo, -|-oo) of Eq. ([1]) may 
be shortened to {t — T,„,f), see Appendix. 

A quasi static approximation of the firing statistics 
is possible when the stimulus s(t) evolves slowly in in- 
tervals of size Tin- In this case, inside all the integration 



hI}-''^ = f dtl... r dtk hk{t-ti,...,t-tk)YliU-ty^ 

= / dxi... dxk hk{xi,...,Xk)Y\{-^iy'- (3) 
Jo Jo tJi 

Replacing Eq. ([2]) in ([1]) and using the definition 
the probability -P[i|s] is expanded as a series of terms 
that only depend on time through the instantaneous 
value of the stimulus and its derivatives at time t 

P[t\s] = Po[s{t)] + Pi[s{t),s'{t)] + P2[s{t), s'it), s"{t)] + ..., (4) 

where, explicitly. 



Pr[s,s', 



" = ' {h,---,jk}/ 



1^- "'^ 



(5) 



and the primed sum over the { ji , . . . , } runs through 
all non- negative integers that sum up to r, that is, 
X^iLi ji ~ (^^^ Appendix for a detailed derivation). 

The value of H-^}''"'-''' is only dependent on intrinsic 
neuronal properties (the filters hk), whereas the param- 
eters \i incorporate the effect of the stimulus s(t). 

For example, the first two terms read 



Po{s) = ho + Y.{\o? Hl^-\ (6) 

+ 00 

Pi{s) = E (^o)''-^ Al [Hl^-<^ + + . . . + . (7) 

fc=i 

Notice that Po{s) is the activation curve: the func- 
tion that relates the firing rate of the cell to the strength 
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of a constant stimulus s. In Eq. ([6]), this activation curve 
is written as a Taylor expansion in s. For arbitrary r, 
each Pr depends on the instantaneous value of the stim- 
ulus and its derivatives: Pq depends only on s{t), Pi 
depends linearly on s'{t), P2 combines quadratic terms 
in s'{t) with linear terms on s"{t), and so forth. There- 
fore, if the stimulus varies slowly in intervals of duration 
Tm, the spiking probability /'[t|s] of Eq. ([4]) becomes 
an instantaneous function of time. This contrasts with 
the original functional of Eq. ([1]) , containing the whole 
stimulus history. 



2.1 Application to linear neuron models 

In a linear neuron model, the firing probability reads 



P[t\s] = ho ■ 



Jt-T„-, 



(t-ti) s(ti) dti. 



(8) 



In this case, the only non-vanishing filters are Hq 
and hi, so H^^''"'-"' is zero for all fc > 1. Hence, the 
expansion ([4]) reads 



Po[t\s\ 
Pl[t\s] 



ho + Ao hO 
Ai Hi 
A2 ^2/2 



(9) 



Pfe[i|s] = Afe H^/k\ 



A linear model is only applicable within a fairly 
narrow range of stimulus fluctuations, given that the 
firing probability must always remain linear and non- 
negative. In Fig. [TJ^ we show two typical activation 
curves, or tuning curves. On the left we see a bell- 
shaped tuning curve, characteristic, for example, of di- 
rection and orientation-selective cells in the middle tem- 



poral visual cortex (|Albrightlll984D . or of motoneurons, 
selective to the direction of arm-reaching movements 



(jGeorgopoulos et al.lll982| ). On the right, the sigmoid- 



shaped function exemplifies the behavior of many rec- 
tifying input-output relations, found for instance in the 
constrast selective cells in the s triate cortex of mammals 
lAlbrecht and Hamiltoiilll982l ). These curves are obvi- 
ously non-linear. In both cases, for small and large stim- 
uli, the firing rate saturates, making subtle stimulus 
discrimination impossible. Hence, a linear approxima- 
tion of the intrinsic neuronal dynamics is only valid for 
stimuli that fluctuate around the shaded regions, where 
the activation curve is well described by a straight line. 
Although linear models cannot be used outside this re- 
stricted range, in many applications they still prove to 
be useful to capture the temporal properties of the cell's 
receptive field. In Fig. [T}3 we show four typical exam- 
ples of filters hi {t) corresponding to ON, OFF, biphasic 



ON and biphasic O FF cells, as reported for example by 
(|Segev et al.ll2006h . or (jCollisch and Meistei][2008al lbl). 
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Fig. 1 Linear neuron models. A: Typical tuning curves. The 
shaded areas represent the range of stimuli where a linear approx- 
imation is valid. B: Typical linear filters hi{t) representing ON, 
OFF, and biphasic (ON and OFF) cells. Spike generation results 
from a convolution of hi{t) with the stimulus. The dotted line 
represents hi = 



In order to remain in the linear portion of the acti- 
vation curve, the parameters Hq and Hi must fulfill 



ho^J2\^'\ 



fc! 



(10) 



where cta^ is the standard deviation of the fc-th deriva- 
tive of the time-dependent stimulus. 

From Eq. ([9]) we see that Po[s(t)] follows the stimu- 
lus linearly, whereas Pi[s{t), s' {t)] introduces a correc- 
tion that is proportional to the derivative of the stimu- 
lus. For the time being, we assume that the higher-order 
terms are negligible. Whenever 3> li/^ |/T,n, 



P[t\s] 



ho - 
ho - 



■H°s{t) + Hls'{t) 



4t) + J^'^'(i) 



ho + H'^sit - 5o), 



(11) 



where Sq ~ —Hi /Hi. The equality between the last two 
lines in Eq. (fTT|) is obtained from the first two orders 
of a Taylor expansion of s{t — Sq). Therefore, if the 
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Fig. 2 Filters and PSTH of 
three linear ON and three 
linear OFF model neurons. In 

the PSTH, gray smooth Unes cor- 
respond to Eq. ([Sj, and black ir- 
regular traces are from multitrial 
Poisson simulations. In cells 2 and 
5, the amplitude of the filters is 
smaller than in cells 1 and 4. In 
cells 3 and 6, the center of mass 
of the filters is displaced to the 
right. In ON cells, the PSTH is 
proportional to the delayed stim- 
ulus, whereas in OFF cells it is 
proportional to the delayed neg- 
ative stimulus. The proportional- 
ity constant depends on the total 
area below the filter. Thus, in cells 

2 and 5, the PSTH has a smaller 
standard deviation than in cells 1 
and 4 (compare the vertical bars 
to the right of the PSTH). The de- 
lay with respect to the stimulus is 
equal to the location of the center 
of mass of the filters. Thus, in cells 

3 and 6 the delay is larger than in 
cells 1 and 4 (compare the slope 
of the lines connecting extremes 
in the stimulus and the PSTH. In 
all cases, ho = 100 Hz and the 
PSTH is measured in [ms"^] 
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linear filter hi{t) is such that |i?j'^|/|i7j'|Tm <C 1, then 
the PSTH follows the stimulus linearly, with a delay 

Hi _ Jo""' fei(f) (^t) dt _ fei(f) t dt 

° HO £"^h^it)dt /;-/M(t)dt • ^ > 

The delay Sq is equal to the center of mass of the 
filter hi(t). Causality imposes that Sq be always posi- 
tive. 

If the filter hi is always positive (i.e., the cell re- 
sponds to depolarizing stimuh, as the ON filter in Fig.[Tj3), 
then Hi is positive, and Hi is negative. The firing prob- 
ability is then proportional to the delayed stimulus, as 
exemplified in cells 1-3, in Fig.[2]A and B. The constant 
of proportionality between the PSTH and the stimulus 
is equal to H^, that is, the area below the filter hi (t). As 
the amplitude of hi{t) decreases (compare cell 2 with 
cell 1, in Fig. [5]), the standard deviation of the result- 
ing PSTH also decreases (compare the vertical bars to 
the right of the PSTHs). In turn, the delay i5o is equal 



to the center of mass of hi{t). Therefore, if hi{t) is dis- 
placed to the right (compare cell 3 with cell 1), then the 
time-lag between the stimulus and the PSTH increases 
(compare the slope of the lines connecting correspond- 
ing points in the stimulus and the PSTH). 

If hi is always negative (i.e., the cell responds to hy- 
perpolarizing stimuli, as the OFF cells 4-6 in Fig. [2] C 
and D), then Hi is negative and Hi is positive. Thus, 
the firing probability is proportional to the delayed neg- 
ative stimulus, as shown in Fig. [21 As the amplitude of 
the filter diminishes (compare cell 5 with cell 4), the 
standard deviation in the PSTH decreases accordingly. 
If the filter is delayed (see cell 6), the time-lag between 
the stimulus and the PSTH increases as well. 

It is important to notice that the amplitude of the 
PSTH is only determined by Hi, and the delay Sq only 
depends on the center of mass of the filter. Other mod- 
ifications of hi{t) that preserve its total area will not 
change the proportionality constant between the stim- 
ulus and the PSTH. In a similar way, distortions that 
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Fig. 3 Comparison between several stimuli and the associated PSTH, for an ON and an OFF linear neuron. 

Left/right scales measure the stimulus {gray)/PSTH (black). Slow Gaussian stimulus: cutoff frequency 50 Hz. Fast Gaussian stimulus: 
cutoff frequency 250 Hz. Ornstein-Uhlenbeck process: correlation time 100 ms. Triangular and square stimuli: up- and downstrokes 
are determined by a Poisson process of mean 25 Hz. Both ON and OFF cells follow the stimulus accurately, except at discontinuous 
points, where an integration roundoff of approximately 5 ms is observed. Notice that the x-scales vary from panel to panel. In all cases, 
ho = 100 Hz and the PSTH is measured in [ms^^] 



do not modify the center of mass of the filter, leave the 
time-lag unchanged. 

In Fig. El the PSTH of an ON and an OFF linear 
neuron is compared to the injected stimulus. Five dif- 
ferent types of stimuli are used. For both cells the delay 
6o is equal to 5 ms. For slow enough stimuli the PSTH 
is equal to the delayed stimulus (or negative stimulus) 
with a time-lag of 5 ms (see slow Gaussian in Fig.|3]). As 
the stimulus varies faster (see fast Gaussian in Fig. [3]) , 
the PSTH becomes a smoothed version of the delayed 
stimulus. In the case of the Ornstein-Uhlenbeck process, 
only high frequencies are filtered out. Thus, for the tri- 
angular wave, the sharp edges of the stimulus at each 
switching point are converted into 5 ms rounded corners 
in the PSTH (5 ms is also the total width of the kernel 
hi). Similarly, the discontinuous square wave produces 
a continuous PSTH, where the rise and the decay times 
between the up and down states also last for approxi- 
mately 5 ms. 

When hi has both depolarizing and hyperpolarizing 
phases (see, for example, the biphasic filters in Fig.[Tj3), 



then the condition \Hl\ ^ |iJ||/Tni is less likely to hold. 
In the extreme case of a perfectly symmetric biphasic 
filter vanishes. The fluctuations in the PSTH are 
no longer proportional to the delayed stimulus. In fact, 
whenever > |-ffi|/Tm, 



P[t\s] 



where 



ho ■ 
ho 

ho - 



Hi s'{t)+Hl s"{t) 
Hi 



s'{t) + || s"it) 



-Hls'it^Si), 



(13) 



(14) 



Hi /(J"" hi{t) t dt ■ 

Thus, in this case, the PSTH follows the first deriva- 
tive of the stimulus, with a proportionality constant 
equal to Hi and a time lag given by Si, as exemplified 
in Fig. 21 For the slow Gaussian stimulus, the PSTH is 
proportional to the derivative of the delayed stimulus. 
When the stimulus varies faster, as in the fast Gaus- 
sian stimulus, the PSTH smoothes rapid fluctuations 
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Fig. 4 Comparison between several stimuli and the associated PSTH, for an ON-biphasic and an OFF-biphasic 
linear neuron. Left/right scales measure the stimulus derivative (gray)/PSTH (black). Slow Gaussian stimulus: cutoff frequency 50 
Hz. Fast Gaussian stimulus: cutoff frequency 250 Hz. Ornstein-Uhlenbeck process: correlation time 100 ms. Triangular and square 
stimuli: up- and downstrokes are determined by a Poisson process of mean 25 Hz. Both the biphasic-ON and the biphasic-OFF cells 
follow the derivative of the stimulus accurately, except at discontinuous points. The Ornstein-Uhlenbeck process has a discontinuous 
derivative, so the quasi-static approach is not applicable. The derivative of a triangular wave is a collection of constant segments. 
The derivative of a square wave is a sequence of delta-like pulses. Notice that the x-scale varies from panel to panel. In all cases, 
ho = 100 Hz and the PSTH is measured in [ms~^] 



out. In the limit case of the Ornstein-Uhlenbeck pro- 
cess, the derivative of the stimulus is discontinuous, and 
the quasi-static approximation is not applicable. There- 
fore, the black curve in Fig.|3]deviates significantly from 
the derivative of the stimulus (in grey: white noise real- 
ization). For triangular and square stimuli, the PSTH 
detects jumps in the derivative. In the case of square 
waves, the derivative of the stimulus is a sequence of 
delta-functions. Although the quasi-static approxima- 
tion is not applicable in this case either. Fig. |4] shows 
that the model captures the transition times correctly. 



2.2 Application to lincar-nonlincar neuron models 

Consider now a neuron model where the probability to 
generate a spike at time t is 



P[t\s] 



ho ■ 



f 

Jt-T„ 



hi{t - ti) s(ii) dti 



(15) 



where g is a static, 
els have been wideH 



eis nave oeen wiaeiy 
Schwartz et al.l200a 




nonlinear funct ion. These mod- 
used in vision (Rust et al. 20051: 



Brenner et al."2000':'Gollisch and Meistei 



and audition ( Nagel and Doupe.2006. : .Lesica and Grothel 
Equation (fTS]) is not expressed as a Volterra se- 



ries, because the static nature of g allows us to work 
with a compact expression, with no need to expand in 
powers of the stimulus. With this model, when consid- 
ering an ON (OFF) ceU, if Eq. HU]) holds, 



P[t\s] 



fc! 



fc=0 

g[ho + Hf sit -So)] 
g{ho)+g'{ho) Hi s{t - So) + 



+ \ 9" {ho) [Hi s(t-5o)y 



(16) 



If we only consider up to the linear term, the simi- 
larity between Eqs. ([TB)) and ^ implies that linear-non 
linear neuron models are formally analogous to purely 
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Fig. 5 The PSTH of linear-nonlinear neuron models. A slow Gaussian stimulus of 50 Hz cutoff frequency and 10"'^ standard 
deviation is presented to a linear-nonlinear model with the same ON filter as cell 1 in Fig.^{H° = 2.506). Left: Nonlinear function 
g{ho). Right: g'{ho). Middle: PSTH obtained from numerical simulations (black line) and from Eg. 1151 (white line). The three cases a, 
b, and c correspond to ho = 250, 300 and 400 Hz, respectively. The mean PSTH is determined by g{ho), as indicated by the dotted 
lines connecting the left and middle panels. Its SD (represented by the height of the grey areas in the middle panel) depends on g'{ho), 
as indicated by the dotted lines connecting the middle and the right panels. The delay with respect to the stimulus is independent of 
ho and g 



linear models: in ON or OFF cells, the PSTH follows 
the stimulus linearly with a fixed time lag. An analo- 
gous argument can be developed for symmetric biphasic 
neurons, showing that the PSTH follows the stimulus' 
derivative. The only difference with the purely linear 
case relies on the constant terms g{ho) and g'{ha): now 
the proportionality constant between the PSTH and 
the stimulus (or its derivative) depends on the opera- 
tion level /iQ. The non- linear function g, nevertheless, 
has no effect on the time lags 5q and i5i , since the latter 
depend only of the ratios Hl/H^ and Hi/ Hi. 

In Fig. [5] we see the PSTH of a linear-nonlinear neu- 
ron model (central panel), for three different operation 
levels hf). The three values can either correspond to 
three different cells with different spontaneous rates, or 
to one same cell driven with three stimuli with differ- 
ent mean values. As is near saturation, the average 
PSTH approaches its maximum possible value. In ad- 
dition, the derivative g'{ho) — > 0, therefore the variance 
of the PSTH decreases. 

We now discuss the third term in Eq. (fT6|) . This 
contribution is not present in the purely linear neu- 
ron model, and is proportional to the squared stimulus, 
that is, the energy of the input signal. Its effect be- 
comes significant for fairly large stimulus amplitudes. 
The quadratic term may or may not have the same 



sign as the linear term, depending on several factors. 
The scheme in Fig. [6] summarizes the combined effects 
of stimulus fluctuations (upward or downward) and the 
concavity of g at the operation level ho, both for the 
linear and quadratic terms of Eq. (fTH)) . For example, 
if g"{ho) > 0, an ON cell transforms positive stimulus 
fluctuations in enhanced peaks in the PSTH. In this 
case, both the linear and the quadratic terms contribute 
with the same sign. Instead, negative fluctuations ap- 
pear less pronounced, since the linear and quadratic 
terms have opposite signs. As a consequence, the PSTH 
is no longer symmetric with respect to its mean value: 
upward fluctuations are amplified, and downward fluc- 
tuations appear as less pronounced. Instead, if ho is 
such that g" {ho) < 0, positive stimulus fluctuations are 
flattened, whereas negative ones appear as enhanced 
troughs. These effects deviate from the purely linear 
prediction, as exemplified in Fig. [71 In case (a), the 
operation point ho is exactly at the inflection point 
of g{ho). Therefore, the third term in Eq. (|16|) van- 
ishes, and the fluctuations in the stimulus are repre- 
sented symmetrically in the PSTH. However, if ho is 
such that g{ho) has negative concavity (case (b)), neg- 
ative fluctuations are amplified, whereas positive ones 
are flattened (see Fig.jB]). This asymmetry is clearly vis- 
ible from the different heights of the two striped grey 
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Fig. 6 The effect of the con- 
cavity of the activation curve. 

Upward or downward stimulus 
fluctuations may or may not be 
intensifled in the PSTH, depend- 
ing on the concavity of the non- 
hnear function g, at the opera- 
tion level ho. The symbols © and 
Q indicate the contribution of the 
linear and the quadratic terms in 
Eq. III6I I. When both contribu- 
tions have the same sign, the fea- 
ture is enhanced 
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areas representing the amplitude of positive and nega- 
tive fluctuations, in the top PSTH of the middle panel. 
The opposite effect (amplified upward fluctuations and 
flattened downward fluctuations) is observed in case (c). 



3 The inter-spike interval distribution 

In Poisson neuron models, the inter-spike interval (ISI) 
distribution /(r) can be written in terms of the PSTH 
as 



fir\s) = ^J^ " P[t\s] Plt+T\s] e^p^^^l'^^ P[t'\s] dt'j , 

(17) 

where T is the total length of the recording, and N 
is the total number of spikes, N — P[t\s] dt. The 
factors P[t\s] and P[t + t|s] represent the probability 
of spike generation at times t and t + t, respectively. 
The exponential factor is the probability that no spike is 
fired inside the interval [t, t + r). It is important to keep 
in mind that Eq. (jl7p is only valid within the framework 
of Poisson neurons, that is, cell models in which the fir- 
ing probability only depends on the incoming stimulus. 
In particular, all the effects of previous cellular activ- 
ity (refractoriness, adaptation) are neglected. Strictly 
speaking, dynamical neuron models cannot be reduced 
to a Poisson formalism, although in some cases an ap- 
proximate reduction is justified (lAviel and Gerstnei]|2006 ) 
Hence, Eq. ((TT)) should be understood as a property of 



theoretical Poisson models, that is only expected to de- 
scribe the behavior of real neurons for intervals r that 
are large compared to the cellular processes that regu- 
late after-spike hyperpolarization and repolarization. 

If is the ISI distribution of a constant stimulus 
s{t) = .s, then 



Po{s) exp[-r Po{s)] 



(18) 



which, in the limit T ^ 00 tends to the classical ex- 
pression for Poisson spike trains /*(t). Hereafter, we 
assume that T ^ r, so the superior integration limit in 
Eq. dni) is replaced by T. 

Decomposing the PSTH using Eq. (|4]) and keeping 
only the lowest order term, the ISI distribution reads 



/(r|s) 1 dt Po[s{t)] Po[s{t+T)] exp |- y*^" PiMt')] dt'j 

(19) 

To further explore Eq. (fTO]) . we separately consider 
three different regimes, depending on how r compares 
to the memory time constant of the system, Tm. 



3.1 The ISI distribution for small t 

When the ISI distribution is evaluated at a time r that 
is much smaller than the typical time scale of stimulus 
fluctuations, the inner integral of Eq. (fT9)) may be re- 
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Fig. 7 Deviations in the PSTH of linear-nonlinear neuron models, for stimuli with large SD. A slow Gaussian stimulus 
of 50 Hz cutoff frequency is presented (see bottom panel). The SD of the stimulus is SxlO""^. The linear-nonlinear neuron model 
has the same ON filter as cell 1 in Fig. [2] (J^J = 2.506). Left panel; Nonlinear function g{ho). Three different operation points ho 
are selected. The corresponding values of g{ho) set the spontaneous firing rate, as shown by the thin connecting lines joining the 
left and middle panels. In case (a), g{ho) is at its inflection point, so the concavity g"(/io) = 0. In (b), g"(/io) < 0, whereas in 
(c), g"{ho) > 0. The resulting PSTHs are displayed in the middle panel, line conventions as in Fig. [5] In case (a), the upward and 
downward fluctuations in the stimulus are symmetrically represented in the PSTH. Instead, positive fluctuations are flattened in 
(b), whereas negative fluctuations are amplifled. The opposite efi'ect is observed in (c). The differential ampliflcation of positive and 
negative fluctuations can be predicted from the value of the second derivative of g, as shown by the thin connecting lines to the right 
panel 



placed by exp{— TPo[s{t)]}. In addition, P[s{t + T)] 
P[s{t)]. With these approximations, 



/{r|s)« l^^di Po[s(i)] rW(r). 



(20) 



The ISI distribution becomes a temporal average of 
the constant-stimulus distribution /'^ (t) , with a weight 
function that is proportional to the PSTH. Moreover, 
for stationary stimuli, the temporal integral of Eq. (|20p 
can be transformed into an integration in stimuli. In 
Fig. [5] we show how to construct the stimulus density 
p{s) from a time-dependent signal s = g{t), by counting 
the fraction of stimuli that fall in the interval (s, s-|-ds) 



P(«) 



Jo 



P{s,t) dt : 



■g{t)] dt 



TJo 



5[t-g-\s)] 



dt 



E 



(21) 



The derivative dg/dt transforms the temporal dif- 
ferentials into stimulus differentials, i.e., accounts for 
the variable widths of the shaded temporal intervals in 



Fig. [5] The points ti are all the times for which g{t) — s, 
and are located at the center of the vertical shaded ar- 
eas of Fig. [5] The times for which the stimulus varies 
slowly, hence, have a large contribution to p(s). There- 
fore, in the stimulus domain, Eq. (|20p becomes 



firls) 



hi 



ds 



E 



Po(s) r{r) 



j P(S) Po{s) r(T) ds. 



(22) 



Here, the first equality results from changing the 
integration variable from t to s, and collecting together 
all the time points ti that correspond to one particular 
s. The second equality derives from the definition of 
p{s) in Eq. ((2T|) . Equation ([22]) is useful to calculate 
the behavior of /(t|s) for small r, without requiring any 
knowledge of the temporal evolution of the PSTH. Only 
the activation curve Po{s) and the stimulus distribution 
p(s) are needed. In particular, two different stimuli with 
the same distribution p{s) but very different temporal 
characteristics (an example is shown below) lead to the 
same /(r|s). 

In order to obtain systematic corrections to Eq. ([^ 
when r becomes comparable to the time scale of stimu- 
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Fig. 8 Construction of the distribution of stimuli p{s) 
from the temporal evolution s{t). One option is to integrate 
the signal g{t) in time. Another possibihty, is to collapse the signal 
on the vertical axis, giving a weight to the (s, s + ds) proportional 
to the total time for which the signal g{t) falls inside this range. 
If ds is small, this time is a sum of terms |ds/dt|~^, one such 
term for each point where g{t) = s 



lus fluctuations, the functions ^'[i|s] in Eq. (|17p should 
be replaced by several terms of their Taylor expansion 
(not just the first). 



In Fig. [9] A, we show a square stimulus of amplitude 
A and period T. This is a binary stimulus, whose dis- 
tribution p{s) — [S{s — A) + S{s + A)]/2 is depicted on 
the left. Integrating Eq. ([2^ . we obtain 

■^(^) = ^ [■^+ ("^ ■^+) + ''^P ■^-)] ' (24) 

where H± ^ ho±H^ A. In Fig. [9] A, we see that the ISI 
distribution corresponding to a slow square- wave stim- 
ulus has a higher peak and a longer tail than the purely 
exponential distribution of decay rate ho obtained for 
a constant stimulus of the same mean. This is because 
there are two exponentials in Eq. (|24p . One of them 
has a slower decay rate than ho, and is responsible for 
the long tail of the distribution, at large r. The other 
exponential decays faster than ho, and dominates for 
small T. Notice that although Eq. was derived in 
the limit of short t, the red curve reproduces the nu- 
merical data up to r values that are 6 times larger than 
the mean ISI, in this case, 10 msec. 



3.2 The ISI distribution for large t 

If the ISI distribution is evaluated for an interval t that 
is much longer than the typical time scale of stimulus 
variations, then the inner integral of Eq. (|17p contains 
many stimulus fluctuations. The number of spikes gen- 
erated in the interval (t, t + t) can thus be assumed to 
be equal to N t/T, and the exponential factor may be 
replaced by exp(— t N/T). Hence, 

/(t|s) ^ 1 exp (-^ J\t Po[s{t)] PoHt + r)]. (23) 

The ISI distribution is exponential in r, and pro- 
portional to the autocorrelation of the spike train. 



3.3 The ISI distribution for intermediate r 

When T is comparable to the typical time scale of stimu- 
lus fluctuations, the ISI distribution must be calculated 
with the full Eq. ^T7\ or keeping more terms for P[t|s]. 



3.4 Examples for the linear neuron model 
3.4-1 Deterministic input currents 
Periodic square wave 



Periodic triangular wave 

We now consider a triangular stimulus of amplitude 
A and period T, as in Fig. [5]B. In this case, p{s) = 
1/2 A, for s e [-A, A]. By integrating Eq. ([22]), we get 

- {nl + 2H+/T + exp (-T ■H+) ] . (25) 

This ISI distribution also contains two exponentials 
with slightly different decay times, producing a higher 
peak at the origin, and a longer tail at large r values, 
as shown in Fig. [5] B. The effect is less noticeable as 
with the square wave, given that this stimulus spends 
less time in the extreme values. 

Sinusoidal wave 

For a sinusoidal stimulus of Fig. [5] C, the stimulus 
density is p{s) = {7rZ\^l - sVZ\2)-i, for s € [-A,A]. 
The numerical integration of Eq. (|22p is depicted in 
the right panel of Fig. [9] C. Also in this case the ob- 
tained distribution differs from the exponential distri- 
bution, by displaying a higher peak at zero, and a longer 
asymptotic tail. Once again the approximation for small 
r values correctly reproduces the entire range of ISIs. 

3.4-2 Stochastic input currents 
The Gaussian stimulus 
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Fig. 9 Examples of ISI 
distributions obtained 
by a linear ON cell, 
driven by slow deter- 
ministic and stochas- 
tic stimuli. Left pan- 
els: stimulus distributions 
p(s). Middle panels: tem- 
poral evolution of the stim- 
ulus. Right panels: ISI dis- 
tributions obtained from 
Eq. ^ for short r (red 
continuous line) and his- 
tograms obtained by nu- 
merical simulation of the 
dynamical system (stairs 
in black line, with grey- 
shaded area). The dashed 
line serves for comparison 
with the exponential dis- 
tribution, corresponding to 
a constant stimulus of the 
same mean value. Square, 
triangular and sinusoidal 
waves: A = 0.025 [stim], 
T = 200 ms, ho = 100 Hz, 
H° = 2.506 ms-l[stim]-l. 
Gaussian signals (d-e) have 
SD = 0.01 [stim] 
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(b) Triangular wave with deterministic switching points 
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(c) Slow sinusoidal wave 
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(d) Slow Gaussian signal (f^ = 50 Hz) 
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For a Gaussian stimulus of standard deviation a and 
cutoff frequency /c 



1 {{aH\{r-r^)X 



TO 



1} 



exp 



4 ^ 



ro) = 



(26) 



where tq = ho/ {aHi)^ . The ISI distribution for small r 
does not depend on the cutoff frequency of the stimulus, 
that is, on its temporal correlations. These correlations 
only participate in determining the range of validity of 
Eq. ([26|) . In the right panel of Fig.[9]D, however, we see 
that the range of validity of Eq. ([26ll extends through- 
out the whole range of ISIs. Hence, the ISI distribution 
of slow Gaussian signals is independent of the cutoff 
frequency. 

In Eq. (|26|). the ISI distribution has a quadratic de- 
pendence on r (both in the factor and in the exponen- 
tial term) . The conditions of validity of the approxima- 
tion {ho ;2> \aHi\) imply that tq must be large. In this 
limit, Eq. (|26p reduces to a purely exponential func- 
tion, visually indistinguishable from the one obtained 
for constant stimulation, as shown in the right panel of 
Fig.[9]D. 

In the limit of large t, the ISI distribution reads 



fir) 



exp (-r N/T) 
TV 



■{H 



Jo 



s{t) s(t + t) dt 



(27) 



The integral in Eq. (P7)) is the autocorrelation of 
the stimulus. For white Gaussian stimuli with cutoff 
frequency /c, the integral can be calculated analytically, 
and 



fir) 



exp {-ho t) 
ho 



hi ■ 



(jj0\2 ^ sin(27r/c r) 



(28) 



where we have replaced N/T by ho- 

The tail of the ISI distribution oscillates with the 
cutoff frequency of the stimulus. These oscillations, how- 
ever, are masked by the rapid decay imposed by the ex- 
ponential term, and the r in the denominator of Eq. (|28[) . 
They are therefore invisible in Fig.[9]D. 

The Ornstein-Uhlenbeck process 

Consider an Ornstein-Uhlenbeck process with a cor- 
relation time A defined by the equation 



A ds/At = -s{t) + cr ^(t), 



(29) 



where ^{t) represents a white noise realization of zero 
mean and unit variance. The distribution of the stim- 
ulus is a Gaussian function of zero mean and standard 



deviation a. Therefore, the ISI distribution for small r 
coincides with Eq. (pS)) . 

In order to obtain the ISI distribution for large r, 
we write Eq. as 



/(,)^_L exp(-^ 



hl + {H^f^ f dts{t) s{t + T) 
^ Jo 

exp(-ftor) r;,2_^(^0)2^2^-./A/2]. (30) 



ho 

As in the previous case, this distribution is domi- 
nated by the exponential term (see Fig. [HlE). 



4 Conclusions 

We derived analytical expressions for the PSTH and 
the ISI distribution of arbitrary Poisson neurons, when 
driven with slow stimuli. Our results depend on a com- 
bination of intrinsic neuronal properties (captured by 
the Hf} '"'•"'), and the external stimulus (represented by 
Xj). In the particular case of linear neuron models, the 
PSTHs of ON and OFF cells are proportional to the 
delayed stimulus. The constant of proportionality de- 
pends on the total area under the linear filter, whereas 
the delay is solely determined by its center of mass. For 
biphasic linear neurons, the PSTH is proportional to 
the delayed stimulus derivative. The constant of pro- 
portionality and the delay can be easily determined by 
performing an elementary integration of functions in- 
volving the linear filter. 

When driven with small-amplitude stimuli, LNL neu- 
ron models can be linearized around their spontaneous 
rate. The resulting PSTH is formally equivalent to the 
one that would be obtained from an effective linear 
model. For larger stimulus amplitudes, however, up- 
ward and downward stimulus defiections are not sym- 
metrically represented in the PSTH. This asymmetry 
is the footprint of the influence of the curvature of the 
activation curve. Depending on whether the activation 
curve operates below or above its inflection point, pos- 
itive or negative stimulus curvatures are selectively fa- 
vored. 

These results might be useful in several applications. 
They could be useful, for example, in the analysis of the 
information obtained by a population of LNL neurons. 
Several studies have aimed to determine the optimiza- 
tion principle that underlies the w ay in which retinal re- 
ceptive fields tile the visual space ( De Vrics and B avloi 

1997 tlVincent et al.ll20oi:ISegev et al.ll2006. : Borghuis et al 



2008HGauthier et al.ll2009t ). These studies focus on the 
population distribution of the spatial properties of re- 
ceptive fields. In this context, it would be interesting 
to understand the population distribution of also the 
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temporal properties of receptive fields. These proper- 
ties involve specific time scales (slow, medium and fast) 
where different time-scales are represented with differ- 
ent frequencies in ON, OFF, biphasic ON and bipha- 



sic OFF cells (|Segev et al.ll200a iGollisch and Meistei 



2008al) . There is presumably also an optimization prin- 
ciple that has guided the evolution towards the par- 
ticular distribution of temporal receptive field proper- 
ties found in each species, probably dependent on the 
statistics of natural images. Our study provides a sim- 
ple picture of how the PSTH of a LNL model depends 
on the incoming signals. Therefore, it should be useful 
to determine the type and amount of information that 
is encoded by a population of neurons characterized by 
a specific distribution of temporal receptive field prop- 
erties. 

In the second part of the paper, we derive the depen- 
dence of the ISI distribution on both neuronal and stim- 
ulus properties. When a Poisson cell is driven with a 
constant stimulus, an exponential distribution is found. 
Here we showed that for variable stimuli, the ISI dis- 
tribution can be obtained by weighting the exponential 
distribution with the stimulus density. Strictly speak- 
ing, this result is only valid in the limit of small ISI. 
However, numerical evaluation of the ISI distribution 
for the case of slow stimuli showed that this expression 
reproduces the whole range of ISIs. Under this approxi- 
mation, the ISI distribution is solely determined by the 
stimulus distribution p(s), independently of the stim- 
ulus' temporal properties. In particular, the ISI dis- 
tribution obtained from a periodic stimulus does not 
depend on the stimulus frequency. In the same way, 
all Gaussian stimuli give rise to the same ISI distri- 
bution, no matter the cutoff frequency. Indeed, Gaus- 
sian stimuli and Ornstein-Uhlenbeck processes produce 
also the same distribution. The obtained densities show 
deviations from the exponential behavior. Specifically, 
the probability density of both small and large ISIs 
is incremented, at the expense of the probability den- 
sity of intermediate ISIs. This result was more visi- 
ble for deterministic stimuli than for stochastic signals. 
Within the class of deterministic stimuli, the largest 
effect was observed in the case of the square wave. Pre- 
sumably, those input currents whose distribution p{s) 
differs most markedly from the constant case (where 
p{s) — S{s — So), for a given sq) are the ones whose ISI 
distribution is most strikingly different from the expo- 
nential function. 

In ( Urdapilleta and Samengol 2009^ . we derived the 
ISI distribution of a perfect integrate-and-fire model 
neuron, when driven with slow stimuli. The analytical 
expression coincided with Eq. (PO)) . implying that the 
result that we obtained here for the particular case of 



a Poisson neuron is also applicable to other dynamical 
neuron models. 

We have focused on slow driving signals. As men- 
tioned in the introduction, there is a myriad of mod- 
ulating processes in the brain that operate in a slow 
time scale. In sensory systems, the biological relevance 
of extremely fast stimuli is questionable, since it has 
been argued ( Atick 19921 ) that the temporal spam of 
receptive fields has evolved as to average out fast stim- 
ulus fluctuations, often contaminated with high levels 
of noise. However, a description of the effect of stim- 
uli whose typical time scales are comparable to those 
of the receptive fields would be highly desirable. This 
would mean to extend the present results to the case 
where the stimulus fluctuations are comparable to t^. 
Presumably these time scales are the most relevant for 
accurate behavioral performance. 
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5 Appendix 

Starting from Eq. ([1]) and replacing (— cx),+oo) integration limits by {t — r„i,i) we obtain 
P[t\s]^ho+ f dti hi{t - h) sih) + f dti / dt2 h2it-h,t-t2) sih) s{t2) + --- + 

t pt 

Ah... dtkhk{t^h,...,t~tk) s{ti)...s{tk) + ... (31) 



Replacing the Taylor expansion of the stimulus, Eq. ([2]), in Eq. ((3T|) . 

+ 00 fc 

P[t\s]=hQ + J2 dti... dtk hk{t-ti,...,t-tk)Y[ s{U) 

+00 „t „t k +00 ^ 

= hQ + Y, Ah... dtkhk{t-h,---,t^tk)\{^^{h~ty' 

= ho+Y, dh... dtkhk{t~h,---,t-h) J2 I[^(ti-*y'' (32) 

where the last sum runs over all the sets of k non- negative numbers ji , . . . , . We now rearrange the sum in 
{ji, . . . ,jk} in increasing order of the product of derivatives. To that end, we group together the terms with the 
same value of the sum r — Jii ^^lat is. 



p[t\s] = ho+Y, dh--- dtk hkit-h,---,t-tk)Y, E li^iti-ty' 

= '^a + J2Y. E n-T/ dti.../ dtkhkit-h,...,t-tk)l[iU~ty 



+ 00 +00 / \ 



= ^o+EE E nf (33) 

where the primed sum j^} '^^'^s over non-negative integers {ji, . . . ,jk} that sum up to r, and Hj}''"'-"' was 

defined in Eq. ([3]). The different terms in Eq. ^ correspond to the different r- values in Eq. 
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